%
% Table 2, Panels A and B. Variance Decompositions
%

clear all
clc

warning('Off','all') ;

path(pathdef) ;

addpath('./routines')
addpath('./routines/solmat')

%% Data

load mhall_01-Apr-2021

SET.EST.burn_in = 1/2;

params_x_ = [] ;
 
for runs=2:maxproc
    params_x_ = [params_x_ ; params_x(round(length(params_x)*SET.EST.burn_in):end,:,runs)] ; %#ok<AGROW>
end

for ii=1:length(params_x_(1,:))
    [a,b]=ksdensity(params_x_(:,ii),'NumPoints',200);
    est_params(ii) = b(a==max(a)) ;
end

save est_params est_params

clear
clc

setup_model

%% Estimation parameters

trend_dates = 1940:0.25:2070 ; 

%convert to quarterly
year_vec  = 1940:2142 ;
datevec_q = 1940:0.25:2142 ;

start_date = 1986 ; SET.start_date=start_date;
end_date   = 2019.75; SET.end_date=end_date;

trend_start = 1960 ;
trend_end   = 2060 ;

trend_dates = trend_dates(find(trend_dates==trend_start):find(trend_dates==trend_end)) ;
SET.trend_dates = trend_dates ;

SET.usecore = 0 ; 
SET.unant_demographics = 0 ;

load_data
data = Zobs;

SET.ss   = length(data(1,:)) ;
nobs     = length(data(:,1)) ;
SET.nobs = nobs ;

SET.EST.hidx      = 1:SET.nobs ;
SET.EST.hidx_no_i = 2:SET.nobs ;

SET.EST.obs_ = [
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve] ;

load ../olg/adj_factors.mat

phi_adj =  phi_s/1 ;
v_adj   =  v_s/1 ;
b_adj   =  A_s/1 ;
a_adj   =  B_s/1 ;
lgv_adj = lgv_vec ;

lgv_adj = [lgv_adj lgv_adj(end)*ones(1,length(year_vec)-length(lgv_adj))];
phi_adj = interp1(year_vec,(1+phi_adj).^1-1,datevec_q,'linear','extrap') ;
v_adj   = interp1(year_vec,(1+v_adj).^1-1,datevec_q,'linear','extrap') ;
b_adj   = interp1(year_vec,(1+b_adj).^1-1,datevec_q,'linear','extrap') ;
a_adj   = interp1(year_vec,(1+a_adj).^1-1,datevec_q,'linear','extrap') ;
tauv    = interp1(year_vec,tauv,datevec_q,'linear','extrap') ;
lgv_adj = interp1(year_vec,(1+lgv_adj).^1-1,datevec_q,'linear','extrap') ;

% restrict changes to a date range
phi_vec = phi_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
v_vec = v_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
b_vec = b_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
a_vec = a_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
tau_vec = tauv(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;
lgv_vec = lgv_adj(...
    find(datevec_q==trend_start):...
    find(datevec_q==trend_end)) ;

trend_params = ...
    [phi_vec(1) ; v_vec(1) ; b_vec(1) ; a_vec(1) ; tau_vec(1) ; lgv_vec(1)] ;

save trend_params trend_params

% store vectors for construction of likelihood
SET.tv_str_chg.ant_chg.phi_vec = phi_vec ;
SET.tv_str_chg.ant_chg.v_vec   = v_vec ;
SET.tv_str_chg.ant_chg.b_vec   = b_vec ;
SET.tv_str_chg.ant_chg.a_vec   = a_vec ;
SET.tv_str_chg.ant_chg.tau_vec = tau_vec ;
SET.tv_str_chg.ant_chg.lgv_vec = lgv_vec ;

SET.tv_str_chg.flag = 1 ;

Tbstar = length(SET.tv_str_chg.ant_chg.phi_vec) ;
SET.tv_str_chg.Tbstar = Tbstar ;

SET.tv_str_chg.ant_chg_params   = ...
    {'phi','v','a','b','tau','lgv'} ;
SET.tv_str_chg.unant_chg_params = ...
    {} ; 

gen_set_params(SET,dyn_in) ;

SET.tv_str_chg.cur_t = 1 ;

TT.t_f = zeros(1,length(data)) ;

mats = tvmats(SET,TT) ;

Qt = mats.Qt ; Qt_trendonly = Qt ;
Gt = mats.Gt ; Gt_trendonly = Gt ;

y_ss = (eye(n_-1) - Qt(1:n_-1,1:n_-1,2)) \ Qt(1:n_-1,end,2)  ;
y_ss = [y_ss ; 1] ; % Constant

x_t = y_ss ;
    
for tt=2:Tbstar
    x_t(:,tt) = Qt(:,:,tt) * x_t(:,tt-1) ;%+ Gt(:,:,tt)*n_T(:,tt) ;
end

model_trend = x_t(SET.EST.obs_, 1:end-2) ;

%%

TT.t_f = T_f ;

mats = tvmats(SET,TT) ;

Qt = mats.Qt ; Qt_trendonly = Qt ;
Gt = mats.Gt ; Gt_trendonly = Gt ;

SET.EST.nparam_est = 11 ; 

SET.tv_str_chg.Qf           = Qt(:,:,SET.trend_dates==(SET.end_date+.25)) ;
SET.tv_str_chg.str_mats     = mats.tv_str_chg.str_mats(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
SET.tv_str_chg.mat_i_f_zlb  = mats.tv_str_chg.mat_i_f_zlb(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
SET.tv_str_chg.zlb          = mats.tv_str_chg.zlb(:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date)) ;
    
Qt = Qt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));
Gt = Gt(:,:,find(SET.trend_dates==SET.start_date):find(SET.trend_dates==SET.end_date));

V  = zeros(nobs,nobs) ; 
H  = zeros(nobs,n_-1) ; 

Ct = Qt(1:n_-1,end,:) ;
Pt = Qt(1:n_-1,1:n_-1,:) ;
Dt = Gt(1:n_-1,:,:) ;

for ii=1:nobs ; H(ii,SET.EST.obs(ii)) = 1 ; end

[x_T, n_T] = ks_(SET, Ct, Pt, Dt, H, V, data, TT) ;

SET.tv_str_chg.cur_t  = 1 ; 
TT.t_f                = 0*T_f ;

mats = tvmats(SET,TT) ;

Qt_nozlb = mats.Qt ; 
Gt_nozlb = mats.Gt ; 
    
Qt_nozlb = Qt_nozlb(:,:,find(SET.trend_dates==SET.start_date):(SET.maxchk+find(SET.trend_dates==SET.end_date)));
Gt_nozlb = Gt_nozlb(:,:,find(SET.trend_dates==SET.start_date):(SET.maxchk+find(SET.trend_dates==SET.end_date)));

%% Table for Variance Decomposition

clc

table=[] ;

for ii=1:6
    table = [table 100*oo_.conditional_variance_decomposition(:,1,ii)];
end

table = table([...
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve
    SET.variable.y
    SET.variable.c
    SET.variable.inve],:) ;

disp('') ;
disp('Variance Decomposition, 2Q horizon') ;
latex(table, '%.0f', 'nomath' ) ;

table = [] ;

for ii=1:6
    table = [table 100*oo_.gamma_y{7}(:,ii)];
end

table = table([...
    SET.variable.iobs
    SET.variable.pinfobs
    SET.variable.dy
    SET.variable.dc
    SET.variable.dinve
    SET.variable.y
    SET.variable.c
    SET.variable.inve],:) ;

disp('') ;
disp('Variance Decomposition, Unconditional') ;
latex(table, '%.0f', 'nomath' ) ;
